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We report extensive simulations of the relaxation dynamics of a self-avoiding polymer confined 
inside a cylindrical pore. In particular, we concentrate on examining how confinement influences 
the scaling behavior of the global relaxation time of the chain, r, with the chain length N and pore 
diameter D. An earlier scaling analysis based on the de Gennes blob picture led to r ~ TV 2 !) 1 / 3 . 
Our numerical effort that combines molecular dynamics and Monte Carlo simulations, however, 
consistently produces different r-results for TV up to 2000. We argue that the previous scaling 
prediction is only asymptotically valid in the limit TV 2> D 5//3 2> 1, which is currently inaccessible 
to computer simulations and, more interestingly, is also difficult to reach in experiments. Our 
results are thus relevant for the interpretation of recent experiments with DNA in nano- and micro- 
channels. 



I. INTRODUCTION 

Polymer chains immersed in solution are subject to 
constant molecular collisions and restlessly undergo con- 
formational changes. Since the motion of each monomer 
or chain segment on a polymer is influenced by the rest, 
the polymer shows unique dynamical properties com- 
pared to those of simple molecules. Confinement in a 
pore or slit can change both the static and dynamic prop- 
erties of a polymer chain qualitatively-^. In the presence 
of cylindrical confinement, for instance, the pore diame- 
ter enters as an additional length scale. An early scaling 
approach* shows how this alters the relaxation dynamics. 

One of the key concepts concerning polymer dynam- 
ics is the rate at which the modes relas^^. Crudely 
speaking, this is a measure of the relative importance of 
each mode in describing the internal motion of a chain 
molecule. In particular, the longest relaxation time or 
the global relaxation time r is of practical importance, 
as it determines how fast the chain reaches its equilib- 
rium. Beyond this time, the internal motion plays no 
significant role, and the chain behaves just like a simple 
molecule. 

A simple model for a polymer is the "ideal" chain, 
which is a phantom chain without self-avoidance. The 
Rouse model provides a description of its relaxation dy- 
namics in an "immobile" solvent without hydrodynamic 
effects, and has been extended to the case of a free self- 
avoiding chain in the presence or absence of hydrody- 
namic effects2^£. In a similar spirit, the scaling behav- 
ior of a self-avoiding chain trapped in a pore of diameter 
D has been studied^; the trapped chain is viewed as a 
linear string of compression blobs^^ of diameter D. 

The problem of a self-avoiding polymer under cylin- 
drical confinement has relevance in a variety of contexts: 



DNA manipulations in nano- or micro-channels^, bac- 
terial chromosome segregation^, and polymer translo- 
cation through narrow pore a 1 ^ 11 , to name a few. Addi- 
tionally, the motion of such a polymer is reminiscent of 
reptation in concentrated polymer solutions^. Conse- 
quently, this problem has been investigated in a number 
of simulation studies£ ' 12 i 13 ' 14 > 15 before. However, to our 
knowledge, the asymptotic blob scaling prediction- for r 
has never been directly confirmed. 

In this work, we perform extensive simulations of a 
cylindrically confined chain, focusing on the slowest re- 
laxation time t. The main goal is to examine dependence 
of t on chain length TV and pore diameter D. For sim- 
plicity, we ignore hydrodynamic effects. Even with this 
simplification, the analysis of r is quite nontrivial, as de- 
tailed below. In an effort to present a concrete picture, 
we combine molecular dynamics (MD) and Monte Carlo 
(MC) simulations. While the MD simulations permit a 
more direct probe into polymer dynamics, the MC sim- 
ulations allow us to consider longer chains. 

Despite the fact that we investigate a wide parame- 
ter space, we observe only an "intermediate" regime in 
our simulations, in which r can be described by power 
laws with a much stronger £>-dependence and a weaker 
TV-dependence than the blob-scaling prediction. Inter- 
estingly, finite-size effects are more sensitively reflected in 
the dynamical quantity r; the static intermediate regime 
is much narrower than the dynamic counterpart. This 
intermediate regime might be relevant for the interpre- 
tation of single-molecule experiments on confined poly- 
mers, e. g., DNA in a nano- or micro-channel^. 

This article is organized as follows. In Sec. HH we 
present a simple physical picture for describing the scal- 
ing behavior of the global relaxation time r of a chain 
trapped in a cylindrical pore. The simulation methods 
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FIG. 1: Cylindrically confined self-avoiding chain. In the 
blob picture, the chain is viewed as a linear string of compres- 
sion blobs depicted as spheres. Inside each blob of diameter 
D, the effect of wall confinement is considered to be minor; 
beyond D, it tends to align the chain in the longitudinal di- 
rection, i. e., along the z axis in our convention. For small 
deformations, the blobs are assumed to be deformed indepen- 
dently of each other, as illustrated by the series of springs in 
the lower figure. 



employed in this work are described in Sec. IHII Sec- 
tion IIVI is devoted to a detailed discussion on the results 
of the MD simulations and their analysis. In Sections 
rVTl and IVII1 we present results of our lattice (MC) sim- 
ulations of the end-to-end distribution of the confined 
chain, from which we extract r and compare it with the 
MD result presented in earlier sections. Finally, we dis- 
cuss in Section IVIIII the implications of our results for 
single-molecule experiments on confined polymers. 



II. THEORETICAL BACKGROUND: 
MODEL AND BEYOND 



BLOB 



Consider a linear self-avoiding chain consisting of N 
monomers of size a, trapped in a cylindrical pore of di- 
ameter D (see Fig. [TJ. In our convention, the z-axis 
coincides with the symmetry axis of the cylinder. If the 
longitudinal chain extension L is much larger than D and 
D 3> a, we can use the de Gennes blob model, which 
views the elongated chain as a string of self-avoiding 
compression blobs of diameter 

Within a blob, the chain is expected to resemble a free 
or unconfined self-avoiding chain of g beads. By equating 
D with the "Flory radius" of the blob^ii! that scales as 
ag u , one obtains g ~ (D/a) 1 ^ , wherein v w 0.59. Linear 
(i. e., longitudinal) ordering beyond D demands that the 
total equilibrium or average chain extension L eq should 
scale as 



L r . 



(L) ~ {N/g)D ~ Na» D 1 



(1) 



where (...) denotes an ensemble average. This scaling 
behavior has recently been confirmed in a theoretical ap- 
proach^. 

To estimate the longest relaxation time r of the chain, 
we assume a Hookian-spring-like response to small exci- 



effective spring constant fc e g of the chain, the relaxation 
time is then reciprocally obtained as r ~ N/k c g, where 
we have assumed that the chain friction is additive. For 
such small excitations, the spring constant of each blob 
should beAi2£ kj, ~ ksT/D 2 ; here and in what follows, 
fee is the Boltzmann constant and T is the temperature. 
We thus obtain 



k R T 



bs 



(2) 



where iVbiobs = L cq /D — N / g is the total number of 
blobs in the chain. 

The mapping of the confined chain onto a series of 
Hookian springs used in Eq.[2]can be understood in anal- 
ogy with a Rouse chain - beyond D, the confined chain 
resembles a one-dimensional Rouse chain 1 ' 6 ... Slow Rouse 
modes are analogous to a random walk in a harmonic 
potential whose spring constant varies as 1/N. The only 
difference is that in the confined case iVbiobs should be 
used in place of N, because now the blob-size D is the 
smallest length scale for one-dimensional "Rouse defor- 
mations" . 

Then, the global relaxation time is 
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where £ is again the friction constant of each monomer. 
This derivation recovers the well-known result of 
Brochard and de GennesA based on the blob picture. 
Hereafter, for simplicity, we use the following conven- 
tions: C, — l,a — 1,/cbT = 1. Then Eq. ([T]) becomes 
L cq = ND 1 - 1 '" w ND-°- 7 , while Eq. © reduces to21 
r~N 2 D 2 -i w N 2 D 3 . 

The effective spring constant k c s of a chain (or any 
elastic rod) can be obtained from its force-extension re- 
lation, from which the relaxation time r can be deduced 
(Eq. [3]). Alternatively, one can relate k e g to the equi- 
librium distribution of the end-to-end distance L, de- 
noted by p(L). The free energy cost for a small change 
in the chain length is a long-time potential effectively 
felt by the chain and determines the equilibrium distri- 
bution p(L), which enables us to relate fc e ff to p(L). For 
L ; 



L cq = (L), the distribution p(L) is Gaussian 



P(L) = 



2tt<j l 



exp 
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(4) 



where ctl is the variance of the distribution. This nat- 
urally arises from the assumption of a Hookian-spring 
like response; the free energy should be ~ (L — L cq ) 2 
for L w L cq . For large deformations, i. e., L <C L 



however, the distribution does not have to remain sym- 
metrical with respect to the line L — L cq , as assumed 
in Eq. HJ see Sec. [VI] for details. The free energy cost 
for a global deformation is then —kgThip(L), and one 



tations from the equilibrium length L cq . If we know the can thus establish fc e fr — (<9 2 /<9L 2 ) [/csTlnp(L)] = a L 
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which agrees with the equipartition theorem 2 ^. The re- 
laxation time r can then be related to the static quantity 
<jl as 



N 



Nai 



(5) 



Thus, the computation of r boils down to that of <7£. 



III. SIMULATION TECHNIQUES 

For the direct, dynamical measurements of the relax- 
ation times, we employed an off-lattice Molecular dynam- 
ics (MD) simulation scheme together with a Langevin 
thermostat for a bead-spring model; to determine static 
properties such as the end-to-end distance distribution, 
we performed a lattice Monte-Carlo (MC) simulation 
that allows for efficient sampling of chain conformations. 
Below, we describe both simulation techniques. 



A. MD details 

In the MD simulations, we used a bead-spring model 
of polymers, which were trapped inside a long cylindri- 
cal tube with a circular or square cross section (we re- 
fer to the latter as "brick-shaped" confinement). The 
bead-bead and bead-wall interactions were modeled by 
the Weeks-Chandler- Andersen (WCA) potential 2 ^ (i.e., 
the repulsive part of the Lennard- Jones potential): 



UwCA(r) = £wca 



(6) 



for r < \f2a and otherwise. Here r denotes the center- 
to-center distance between two beads, or the distance of a 
bead center from the confining cylinder minus a (i.e., the 
cylinder wall excludes the centers-of-mass of the beads) . 
With our choice of ewcA = 1 ksT, this potential mod- 
els soft beads of diameter a: their centers cannot come 
much closer to each other than a and are confined to 
be on the inside of the confining wall. In the simula- 
tion, a defines the basic length scale and ewCA the en- 
ergy scale, i.e., we choose a — 1 and cwca = lfcsT; in 
addition, we fix the mass scale so that the mass of each 
bead m = 1. This automatically sets our basic time scale 
as twca = ciy m/ewcA = 1- Following the conventions 
in the theory section, we will henceforth omit the units 
of these quantities. 

The bond between two neighboring beads, which en- 
dures chain connectivity, was modeled by the FENE (fi- 
nite extensible nonlinear elastic) potential^ 



U F {r) = --4e F ln 



'-'5 



(7) 



where r is again the distance of the bead centers and ep 
is the interaction strength. In the present simulations, 



we chose e F — 10, which in combination with the our 
WCA potential results in a typical bond length of 1.027. 

We simulated this system using the simulation package 
ESPResSo^. To integrate the equations of motion, wc 
employed a velocity- Verlet MD integrator with a fixed 
time step of 0.01; the system was kept at a constant tem- 
perature by means of a Langevin thermostat with a fixed 
friction of £ = mr^ A = 1 . The effect of the thermostat 
is not only to keep the temperature constant but also to 
ensure that each monomer moves diffusively, rather than 
ballistic. This is to mimic the effects of solvent viscosity 
in the high-damping limit, in which the inertia term can 
be dropped 2 -. 

Initially, the chain was created as a biased, non-self- 
avoiding random walk confined within the pore. To equi- 
librate the system and remove possible bead-bead over- 
laps, we simulated the system for a few thousand steps 
with a truncated WCA potential. In other words, the 
WCA potential was modified such that the potential is 
linear for distances smaller than a certain cutoff radius. 
We reduced this cutoff gradually, until the potential had 
eventually converged to the full WCA interaction. 

Next, we equilibrated the system further for a total 
equilibration time of T eq , until the end-to-end distance 
and the radius of gyration had converged to their steady 
averages. This was followed by the sampling phase with a 
duration of T to tai, during which we periodically recorded 
configurations for analysis. Between each two recorded 
configurations we waited for a time of T samp i c to reduce 
the statistical dependence of the configurations. For 
the parameters of the performed simulations, see Ap- 
pendix [XJ 



B. MC details 

In an effort to sample a wide parameter space, we com- 
plement our MD simulations with lattice-based MC simu- 
lations. In our MC simulations, the polymers were mod- 
eled as a self-avoiding walk on a cubic lattice between 
hard walls forming the cylinder. In contrast to the MD 
simulations, where we typically used a cylindrical pore 
with a circular cross section, we mostly used one with a 
square cross section of area D 2 in the MC simulations. 
The latter is more compatible with the geometry of the 
cubic lattice (on which the polymer lives). Additionally, 
hard-core repulsions were defined as excluding conforma- 
tions where two monomers occupy the same lattice site. 

For the Monte Carlo moves, we adopted the "worm- 
hole" method of Houdayer— . This algorithm consists 
of a reptation-like motion so that the first bead jumps 
through a "wormholc" to a random position — not to be 
confused with reptation as a physical process in concen- 
trated solutions 2 '-. During the wormhole move, the poly- 
mer therefore consists of two disconnected segments. The 
reptation is continued, until a valid (i.e., connected) poly- 
mer conformation is obtained. This algorithm is known 
to be very efficient for sampling the chain conformations 
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FIG. 2: Equilibrium chain extension L eq (in units of Rf) 
of a cylindrically confined self-avoiding chain as a function of 
the relative confinement D/Rf. The roman numbers denote 
the three regimes: the strong-confinement (I), intermediate 
(II), and free chain regime (III). The diagonal dashed line in 
regimes I and II represents the asymptotic power law L cq ~ 
ND~ a ' 7 , the dashed line in regime III corresponds to the free 
chain limit L oq = y/T/3Re. The data points for N up to 300 
are MD results, while those for N > 500 are MC results. 

of polymers with excluded volume in confinement or in 
dense polymer melts. We only use this wormhole method 
to study static properties, since, by its nature, the dy- 
namics based on this method is completely unphysical. 

Between two chain conformations, we performed 1000 
wormhole steps to assure statistical independence, and 
sampled in total 20,000 conformations per parameter set. 
Note that this algorithm requires no equilibration, be- 
cause of the way it is constructed^. 

IV. THE END-TO-END DISTANCE 

If a chain is confined in an open cylinder with a di- 
ameter much smaller than its natural size, the chain will 
be stretched along the cylinder axis. Fig. [5] shows the 
confinement-induced chain stretching L cq / Rp as a func- 
tion of the dimensionless parameter D/Rp. L cq denotes 
the equilibrium end-to-end distance parallel to the cylin- 
der axis, and Rp the Flory radius, i. e. the unconfined 
equilibrium end-to-end distance. In the figure, three dif- 
ferent regimes are identified: 

I. To the region below D w 0.2i?p we refer as the 
"strong" confinement regime, in which our simu- 
lation results are captured reasonably well by the 
scaling prediction L cq ~ ND 1 - 1 ^ ~ ND o r . The 
data for large N tend to follow this scaling bet- 
ter, because for small chains, the tube diameter 
approaches the size of the beads, in which case a 
"blob" only consists of a few beads. 

II. The range between D « 0.2Rp and D w Rp, we 
denote as the "intermediate" confinement regime. 



FIG. 3: Global relaxation time r g of the radius of gyration R s , 
in units of the relaxation time of a corresponding free chain, 
Tg rcc ) as a function of the reduced diameter D/Rp. The solid 
line marks the asymptotic scaling law of r ~ N 2 D 3 . The 
inset shows a graph r g /N 2 vs D, demonstrating r ~ JV 2 for 
small diameters. 

In this regime, we find in our simulations a non- 
negligible chance for the chain to revert its direction 
along the cylinder axis or at least partially fold back 
one of its ends. Therefore, the average L (or L cq ) 
decreases faster with increasing D than what one 
may naively expect from the blob picture. 

III. In the region above D w Rp, the chain is essentially 
free; we call this therefore the "weak" confinement 
regime. Here, the chain extension is practically in- 
dependent of D, and is given by L = yj\/ 3Rp . 

As we will argue below, the intermediate regime also 
appears in the behavior of the relaxation times t: how- 
ever, it extends to much smaller values of D/Rp. Reach- 
ing the dynamic intermediate regime for r therefore re- 
quires much larger chains than reaching its static coun- 
terpart for L cq . 

Note that, for the same number of monomers, the mag- 
nitude of Rp is different for the lattice (MC) and off- 
lattice (MD) models. The Rp obtained with the MD 
simulations is approximately 30% larger than that found 
in the MC simulations. This difference can be absorbed 
into the effective monomer size a, which is proportional 
to Rp. Because L eq scales as a 17 , this difference is more 
pronounced in the estimates of L cq : for the same chain 
length and diameter, the chains in the MD simulations 
are about 60% — 100% larger compared to the MC sim- 
ulations. 



V. RELAXATION TIME MEASUREMENTS 

We now turn to the global relaxation time r of a con- 
fined polymer. In the case of strong confinement, one 
expects the blob-scaling result r ~ N 2 D 2 ~ 1 / 1 ' w N 2 D o:i 
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to be valid, in parallel with the static scaling law L cq ~ 
LD X ~ X I V . For weak confinement or D w one expects 
the scaling behavior of the relaxation time to cross over 
to that of a free polymer—: r ~ N 1+2u ss N 2 - 2 . However, 
little is known about the relaxation time in the dynamic 
intermediate regime. Below we present the relaxation 
times determined in our simulations. As it turns out, we 
are not able to observe the expected blob-scaling for t, 
even for D <C 0.3i?F- This indicates that the dynamic 
intermediate regime spans a much wider parameter space 
than the static intermediate regime. 

We determine the relaxation time t primarily from the 
slowest exponential decay of the autocorrelation function 

c (t) = ((i?g(t)-(i? g ))CR g (0)-(i?g))) () 

~ (Rl) - (Rz) 2 

where R g — (1/iV) J^i ( z i ~ ( z i)) 2 denotes the longitu- 
dinal component of the radius of gyration. The slowest 
exponential decay r g we obtained from fitting the long- 
time behavior of C g (t) ~ e~*/ Tg . r g turns out to be rather 
insensitive to the exact definition of i? g ; the relaxation 
time tl of the end-to-end distance differs from r g by less 
than 20%, and its overall behavior is the same (r^-data 
not shown). We also determined the relaxation time by 
measuring the relaxation of a stretched chain. These re- 
laxation times agree with the times obtained from the au- 
tocorrelation analysis, as expected on the basis of linear- 
response theory. However, the correlation function of 
the radius of gyration allows us to obtain better statis- 
tics. The global relaxation time r g obtained from our 
MD simulations is shown in Fig. [3] 

For D > 0.5i?F, T g quickly approaches that of a free 
chain. However, it has a maximum around D rs 0.5i?F, 
beyond which the relaxation time decays slowly with in- 
creasing D and eventually becomes D-independent for 
sufficiently large D. We argue that this non-monotonic 
behavior is associated with the reversion of the chain di- 
rection or chain back-folding. This process happens fre- 
quently in our simulations for D > i?p , and its timescale 
increases with decreasing D, as we will show below. Be- 
low D = OARp, however, we do not observe any back- 
folding in our trajectories, so that the observed relaxation 
time is indeed dominated by the fluctuations of the chain 
length. 

For D is i.e., in the static intermediate and strong 
confinement regimes, the relaxation times we find in our 
MD simulations approximately scale as 

tmd ~ N^D 1 - 3 . (9) 

This result has a much stronger /^-dependence and a 
slightly weaker TV-dependence than the scaling result, 
which varies as N 2 D 3 . We refer to this as a dynamic 
intermediate regime, which spans deeply into the static 
strong confinement regime. We will present further anal- 
ysis of this behavior based on blob statistics in Sec. IVIl 
In a small range in the static intermediate regime, we 
observe that the relaxation time increases more steeply 
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FIG. 4: Global relaxation time r z obtained from the autocor- 
relation function of L z , i.e., the longitudinal component of the 
end-to-end vector. The r z is given in units of the relaxation 
time ri ree of a corresponding free polymer. Our estimate of t z 
increases strongly with decreasing D, in contrast to r g plotted 
in Fig. El 

with D than indicated by the result in Eq. (0). However, 
this range quickly decreases with increasing chain length, 
and we will not consider it further. Note that in our 
simulations, we cannot consider diameters smaller than 
D = 2. For narrower tubes, two beads can no longer 
pass each other without paying an energy penalty. This 
artifact significantly reduces the relaxation time, as can 
be seen from the anomalous behavior of t z towards the 
smallest D /Rp values in Fig. [3] (see the data points at 
the left end of each curve). Also, in this case, g w 1, and 
thus the blob picture breaks down. 

To illustrate the significance of chain reorientation or 
back-folding in the weak confinement regime, we have 
plotted in Fig. 2] our results for the global relaxation 
time of the longitudinal end-to-end vector, the differ- 
ence L z = Zhoad — z taii of the z-coordinates of the head 
and tail beads (obviously, L — \L Z \). This relaxation 
time, denoted by r z , is much longer than the L- or R g - 
relaxation times; in our simulations, we were only able to 
track it down to D/Rp w 0.4. In the case of a free poly- 
mer, the relaxation times obtained from the fluctuations 
of its end-to-end vector and end-to-end distance are ex- 
pected to be comparable to each other. In the presence 
of confinement, however, the distribution of L z has two 
minima at L z « ±L oq , separated by a free energy bar- 
rier; these two minima correspond to the two possible 
chain orientations. In this case, chain reorientation, i. e. 
"tunneling" from one minimum to the other, represents 
the slowest mode. It requires a complete back-folding of 
the chain, which becomes unlikely with increasing chain 
length. Since the chain length increases with decreasing 
D, the reorientation timescale increases fast with D. 

As an alternative approach, we have determined the 
relaxation time from direct stretch-release simulations. 
In the linear-response regime, the relaxation times ob- 
tained from the autocorrelation function of the chain 
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FIG. 5: Comparison of two relaxation times: the global 
relaxation time r g of the longitudinal component of the ra- 
dius of gyration, obtained from the autocorrelation function 
(filled symbols), and that from direct stretch-release simula- 
tions (open symbols). The dotted lines are guides to the eye. 



length should agree with the results of simulations that 
measure the relaxation of a slightly stretched chain. In 
order to show this, we carried out stretch-release simu- 
lations as follows. We started our simulations with equi- 
librated chains in confinement, whose end beads were 
initially fixed at a distance L = L cq + 2<tl, where a\ 
is the variance of the chain length L. We then released 
the length constraint and recorded the time evolution of 
L, whose long-time behavior shows an exponential relax- 
ation to its equilibrium value L cq : L(t) — L cq w e~ t ^ TL 
for large t. In Fig. [5l we have plotted the resulting tl , 
together with the global relaxation time r g obtained from 
the autocorrelation function. While the two sets of data 
agree quite well with each other, the autocorrelation anal- 
ysis yields better statistics. This analysis demonstrates 
the reliability of our results for r g . 



FIG. 6: MD results for the distribution of the end-to-end 
distance for N = 200 and various pore sizes. The solid lines 
are Gaussian fits to the graphs. 
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FIG. 7: Ratio of the relaxation time r g to Na\ obtained 
from MD simulations, as a function of the diameter D/Rp. 
As evidenced in the figure, this ratio tends to a D-independent 
constant (~ 0.1) as D decreases. 



VI. LINKING BETWEEN DYNAMICS AND 
STATICS: THE L-DISTRIBUTION 

The Molecular Dynamics simulations, that we have 
presented so far, are limited to chain lengths up to 
N = 300, because the relaxation of longer chains in con- 
finement becomes prohibitively slow. To circumvent this 
difficulty and to compliment the MD results, we carried 
out Monte Carlo (MC) simulations. Note that one can- 
not deduce dynamical information directly from the MC 
wormhole algorithm, since the artificial dynamics of this 
method relaxes the system nonphysically fast. However, 
for strong confinement, it is possible to connect the re- 
laxation time, a dynamic quantity, to a static property, 
namely the distribution of the end-to-end distance, as 
illustrated in Sec. HU 

In Fig. [HI we have plotted our MD results for the dis- 
tribution of the end-to-end distance L, denoted by p(L). 



As shown in the figure, this distribution is essentially 
Gaussian - the Gaussian fit works better for small D, 
i. e., in the strong confinement regime, as also indicated 
in Sec.|TTJ For large D, p(L) is only described well by the 
Gaussian in the region on the right side of the peak. The 
noticeable discrepancy towards the left end of the distri- 
bution can be attributed to partial chain back-folding, 
which tends to widen the distribution and creates an ex- 
ponential tail. 

According to Eq. ([5]), the relaxation time is related to 
a\ via r ~ Na\ in strong confinement. This means that 
in this regime, the D dependence of r g can be deduced 
from ox, since both have the same D-scaling. Fig. [7| 
shows the ratio of the relaxation time r g and the vari- 
ance Na\ of the L-distribution, both obtained from the 
MD simulations. As expected, the ratio in the figure 
tends to a D-independent constant (« 0.1) as D/Rp 
decreases. The non-monotonic D-dependence around 
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FIG. 8: MD results for the end fluctuation <j\ of a confined 
polymer, rescaled by N, as a function of D. The peaks are 
all located at D » 0.5-Rf. The solid line represents the De- 
scaling expected from the blob approach. In the inset, <y\ is 
rescaled by JV ' 75 , and the dashed line corresponds to D 1A , 
in accordance with our best fit. 
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FIG. 9: MC results for the end fluctuation erf, of a confined 
polymer, rescaled by N , as a function of D. The solid line 
represents the D°' 3 -scalhig expected from the blob approach. 
In the inset, a\ is rescaled by TV 0,75 , and the dashed line 
corresponds to D ' 75 , in accordance with our best fit. 

D rs 0.3i?F is mainly due to a significant overshoot in a\ 
for D w 0.5i?F- This overshoot is caused by the widening 
of the L-distribution towards small L due to the onset of 
chain reorientation or back-folding. Having established 
the relation r ~ Na^ , we can focus on a\ and investigate 
its behavior more closely. 

In Fig. [H] and O we represent our results for a\ from 
the MD simulations and the MC lattice simulations, re- 
spectively. The blob picture presented in Sec. [TT1 suggests 
that a 2 L scales as ND - 3 , i.e., a\ is extensive. However, 
both simulations indicate that o\ is not proportional to 
N, and has a much stronger D-dependence. More pre- 
cisely, our data can be fitted well by u\ « J\T°' 75 -D 1-3 for 
the MD simulations, and a\ w iyo.75^,0.75 for the MC 
simulations. The main difference between the two sets of 




FIG. 10: Typical blob size distribution for D = 5 and chain 
lengths iV = 100 and N = 200. A plain Gaussian distribution 
(solid line) fits the data well; for small Lb, the tail is better 
described by an exponential (dashed line). 

simulations is therefore the Z?-dependence; we will argue 
below that this is a lattice/off-lattice effect. While the 
MD result is consistent with r in Eq. [9l the MC result 
leads to the estimate 

tmc ~ N^D - 75 . (10) 

In the next section, we examine this unexpected scal- 
ing more carefully. For this, it proves useful to examine 
single-blob statistics, especially the variance <j\ of the 
single-blob size distribution. This will enable us to inves- 
tigate the N- and /^-dependence of r independently. 

VII. ANALYSIS OF SINGLE BLOBS 

To study our system at the level of a single blob, we 
determine the number of beads per blob g from our sim- 
ulations using the relation L/D — N/g (Eq. (JTJ) ) . We 
then view each sub-chain of beads i + 1, . . . , i + g, con- 
sisting of g consecutive beads, as a realization of a blob. 
For example, we sample the distance \z%-\-i — i- e., 

the longitudinal distance between two "end" beads, as 
the "end-to-end" distance Lh of a blob. In the strong 
confinement limit, the average size of a blob constructed 
this way is indeed Lf, D. 

Fig. [10] shows our MD results for a typical blob size 
distribution p(Lb), which is well represented by a combi- 
nation of two functions: a Gaussian distribution around 
the peak and an exponential one for small L^. In this re- 
spect, the blob size distribution resembles a short chain 
in confinement i?p w D, as indicated in Fig. [BJ In both 
cases, the effect of confinement is marginal and chain 
back-folding is noticeable, leading to the exponential dis- 
tribution for small Lh- Furthermore, using both simula- 
tions (MD and MC), we found that the distribution of L\, 
is practically independent of the chain length if L > 3D, 
i.e., if the chain consists of at least 3 blobs. 
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FIG. 11: MD and MC results for the ratio of the variances 
of the longitudinal end-to-end distance and the blob size, as 
a function of the number of blobs iVbiobs = N/g = L eq /-D. 
The solid line represents the extensiveness a\ = afiVbiobs 
expected. Two types of confined space were used: a cylinder 
and a brick, i. e., a tube with a circular and rectangular cross 
section, respectively. 
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FIG. 12: MD and MC results for the average variance erf of 
the single blob size, as a function of the tube diameter D. As 
in Fig. 1111 we used two types of confined space, a cylinder and 
a brick. The solid line represents the theoretically expected 
-D 2 -scaling. 

Fig. QT] shows crj^/a 2 , the ratio between the variances 
of the end-to-end distance and the blob size. If the fluc- 
tuations of the individual blobs are independent, this 
ratio should equal the number of blobs iVbiobs = N/g, 
i. e., the end-to-end fluctuation is extensive. The data 
was obtained from both MD and MC simulations, and 
for two different shapes of the confined space each: a 
cylinder and a brick, i. e., a tube with a circular and 
square cross section, respectively. The discrepancy be- 
tween the four cases is minor, however, we observe clear 
deviations from the linear relationship. The results in 
Fig. [TT] for up to iVbiobs ~ 30 are, in fact, better de- 
scribed by a 2 L /a 2 ~ A^iots- For lar g e N hXohs , the data 
falls onto a line, however shifted up by about 1 and with 



a slope of about 0.9. The shift is a finite-size effect that 
simply reflects that the size fluctuations of the two end- 
blobs are much larger. The slightly lower slope in turn 
indicates that the fluctuations between the blobs are not 
completely independent. We will discuss the implications 
of this finding below. 

To further proceed with our discussion regarding the 
breakdown of extensiveness of the end fluctuation, we 
first focus on the variances a 2 , as shown in Fig. [T^J We 
only took into account data points with L > 3D, for 
which the chain length dependence of Lj, was negligible. 
The difference between the cylinder and the brick is again 
just an insignificant prefactor. However, while the shape 
of the Lb-distribution is very similar for the MD and 
MC simulations (i. e., approximately Gaussian), there 
are significant differences in the D-scaling of a 2 . The 
MC lattice simulations obey the theoretically expected 
a 2 ~ D 2 law (see Eqn. ([2])) reasonably well, while in the 
MD simulations, a 2 rather follows a 2 ~ D 2 6 . 

The most plausible explanation for this difference is 
through effects near the cylinder surface. The predic- 
tion a 2 ~ D 2 requires that the radial bead distribution 
can be written as a function of p/D , where p is the 
radial coordinate 20 . However, in the off-lattice simula- 
tions, we observe layering effects of the width of about a 
bead diameter at the cylinder surface, which are absent 
on the lattice^. Because this surface layer for D = 6 
still contains about 30% of the total volume, this makes 
a significant contribution. Nevertheless, the radial distri- 
bution converges slowly to the universal distribution, and 
we expect the MD result to cross over to a D 2 -law for suf- 
ficiently large D; unfortunately, this is computationally 
inaccessible. 

The single-blob analysis presented in this section is 
consistent with the unexpected D-dependence of a\ 
(and that of the relaxation time): We observe <r\ w 
(N/g)°- 75 al = N 75 D- 13 af, where a 2 ~ D 2 e in the 
MD simulations and a 2 ~ D 2 in the MC simulations. 
In any case, a\ cx a\. This results in a\ « N°- 75 D 1 - 3 
for the MD simulation and a\ m N°- 75 D°- 7 for the MC 
simulation, which are both in good agreement with our 
numerical findings for a\. Finally, when combined with 
our <jl analysis, the relation ([5]) leads to a relaxation time 
t ~ N 1 D 1 ' 3 , again in good agreement with our earlier 
observations. 



VIII. CONCLUSIONS 

We have performed extensive simulations of a self- 
avoiding polymer under cylindrical confinement. In an 
effort to sample a wide parameter space, we have com- 
plimented our off-lattice MD simulations by lattice MC 
simulations, and considered chain lengths N up to 2000. 
Most notably, both types of simulations have produced 
global relaxation times r that have a much stronger D de- 
pendence than expected from the earlier blob approach 1 . 
While our analysis of the longitudinal chain size L for 
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the case D ;$ 0.2i?p tends to support the picture of lin- 
ear ordering of blobs as assumed in the blob approach, 
our results for r for the same case D ;$ 0.3i?F do not 
show the dynamic blob scaling regime, in which the blob 
scaling of r holds. 

Our MD simulation directly measured the time evo- 
lution of the end-to-end distance from which r was 
obtained; in our MC simulation we extracted r from 
the variance of the end-to-end distance distribution a\ 
through the relation r ~ Na^, which we have demon- 
strated to hold in strong confinement. The absence of 
the dynamic blob-scaling regime in both types of simu- 
lations has been attributed to finite-size effects. To this 
end, we have examined the chain at the level of a single 
blob. In particular, we have examined the D-dependence 
of a 2 and examined its relation with a\ . Even with up to 
-Wbiobs = 30 blobs, o~\ = -/VbiobsOfc does not hold. In other 
words, the assumption of extensiveness of ax, which un- 
derlies the blob picture, is not satisfied in our simulations. 
Moreover, the fluctuations of the individual blobs in our 
bead-spring model also suffer from finite size effects; in 
our MD simulations, we could not reach the expected 
of ~ D 2 scaling due to large surface contributions. This 
indicates that the number of beads per blob also has to 
be significantly larger than 100, corresponding to a diam- 
eter which is at least 15 times larger than the effective 
bead diameter. 

Our main findings can be summarized as 

t w N 175 D~ 13 a 2 . (11) 

When combined with our MD result a 2 « D 26 (obtained 
for 5 < N < 300 and 2 < D < 15), this relation leads 
to r w N^D 1 - 3 . If our MC result a 2 « D 2 (obtained 
for 200 < N < 2000 and 2 < D < 20) is used instead, 
we obtain r « iV 1 - 75 !? - 7 . Both are quite different from 
the asymptotic scaling of r ~ N 2 D ' 3 . Nevertheless, our 
results for r seem to hold for 0.05 < D/Rp < 1, and 
should therefore be experimentally observable. 

Importantly, our simulations show that, for reaching 
the asymptotic scaling limit with a bead-spring model, 
about 10,000 beads are necessary. The aspect ratio L/D, 
which is identical to the number of blobs, has to be at 
least 30. The chain-length requirement for the asymp- 
totic scaling regime N ^> D x l v 1 can then be satis- 
fied, however, obtaining the relaxation of such a system 
by simulations will remain an intractable task for some 
time. For smaller chains and diameters, the blob-scaling 
picture has to be used with due caution. This is similar 
to recent findings for a polymer stretched by an external 
force without confinement^, and for the force-stretching 
relation of a polymer in cylindrical confinement^. 

In fact, we notice that even in experiments this regime 
is hard to reach. For instance, consider DNA molecules 
trapped in a cylindrical pore. Although the Kuhn seg- 
ments of DNA are cylindrical and therefore our bead- 
model is not directly applicable, one can obtain some 
rough estimate on the minimally necessary system di- 
mensions. To reach the classical self-avoiding polymer 



limit inside a blob of size D, this diameter has to be 
significantly larger than the persistence length of DNA, 
which is 50nm. Both our simulations and scaling analy- 
sis^, indicate that a pore diameter of about l^m is neces- 
sary. The extension of the chain in the pore should be at 
least 30 times the pore diameter, which requires a chain 
length of more than a million base pairs. This length is 
significantly larger than the 164,000 base pairs used in a 
recent experiment 7 . 
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APPENDIX A: MD SIMULATION PARAMETERS 

In the following tables, we list the simulation param- 
eters used for the different MD runs, namely the simu- 
lated pore diameters D, the equilibration times T eq and 
the sampling times T to tab as wen as t ne time between 
two successive samples T samp i c . The number of samples 
is hence given by Ttotal/^sample- Times are measured 
in multiples of the microscopic time twca, which corre- 
sponds to 100 MD time steps, i. e. the total number of 
performed time steps is equal to 100T to t a i- 
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